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ABSTRACT 

We study the onset of the Marangoni instability in a liquid layer with two free nearly 
insulating surfaces heated from below. Linear stability analysis yields a condition for the 
emergence of a longwave or a finite wavelength instability from the quiescent equilibrium 
state. Using the method of asymptotic expansions we derive a weakly nonlinear evolution 
equation describing the spatiotemporal behavior of the velocity and temperature fields 
at the onset of the longwave instability. The latter is given by AM = 24, AM being 
the difference between the upper and the lower Marangoni numbers. It is shown that in 
some parametric range one convective cell forms across the layer, while in other parametric 
domains two convective cells emerge between the two free surfaces. 
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I. INTRODUCTION 


A quiescent liquid layer exposed to a thermal gradient can exhibit a sudden onset of 
fluid motion under certain conditions. When the liquid layer has one or two free surface(s) 
and the thermal gradient is perpendicular to the free surface(s), one of the instability 
sources is the surface tension dependency on temperature, also know as the Marangoni ef- 
fect. Since Pearson 1 first theoretically investigated the Marangoni instability phenomenon 
in 1958, a great deal of attention has been given to this subject due to its important im- 
plications in various industrial processes, such as the containerless processing of crystals 2 , 
the welding pool technologies 3 , and fluid dynamics in a low-gravity environment 4 . Various 
different aspects of the Marangoni instability have been studied in the literature 1,5-14 . 
However, the great majority of these analyses deal with liquid layers with a bottom rigid 
surface and a top free surface (hereafter referred as the rigid-free case), while little research 
exists on the Marangoni instability in layers with a top and a bottom free surface (hereafter 
referred as the free-free case). 

This physical problem of the free-free case is of significance for understanding the 
Marangoni instability in multilayer liquid systems 15,16 and liquid sheets held between 
two gas phases 17 . In addition, it would enhance the fundamental understanding of the 
Marangoni phenomena. It should be noted that Rayleigh, in his classical work 18 , solved 
the buoyancy-induced thermal instability in an infinite liquid layer for three different top- 
bottom surfaces, i.e., the rigid-rigid, the rigid-free, and the free-free case. The influence of 
the top and bottom surfaces can be clearly seen in the different critical Rayleigh numbers 
for the onset of convection, a rigid boundary having a stabilizing effect. In addition to this 
effect, another factor affecting the surface tension-induced instability in the free-free case 
is that as the liquid layer is differentially heated, the free surface on the cold end becomes 
thermally destabilized while the free surface on the warm end actually provides thermal 
stabilization, if the case of surface tension decreasing with temperature at both surfaces 
is considered. It is the objective of this work to investigate the Marangoni instability in 
the free-free case, to establish the critical threshold for its onset, and to study the steady 
convection pattern slightly above the onset under certain limiting conditions. 

Georis et al 15 studied the Marangoni instability of a tri-layer in microgravity analyti- 
cally in the linear regime and numerically in the nonlinear one under the assumption that 
the values of the Marangoni numbers at both free nondeformable boundaries are equal. 
They showed that the convection is driven by one destabilized surface while the other one 
is stabilizing. Funada 17 investigated the Marangoni instability in a liquid sheet for both 
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the cases of nondeformable and deformable free surfaces, assuming that the Marangoni and 
Biot numbers are the same at both surfaces. The critical value of the Marangoni number 
was found. In particular, he showed that in the case of nondeformable free boundaries 
the instability is to perturbations of a finite wavelength, small wavenumber (longwave) 
perturbations being stable, and the critical Marangoni number depends on the Biot num- 
ber. It was also found that another type of instability may emerge if the boundaries are 
deformable. 

In this paper we study the more general problem in which the ambient phase can 
be different at the two surfaces. We find that under certain conditions (i.e. for small 
Biot number at both surfaces and for a certain range of ratio of the Marangoni numbers 
at the two surfaces) there is a longwave instability at onset. Therefore, just as in the 
rigid-free case for small Biot number 8,9,11 , asymptotic methods are applied to study the 
weakly nonlinear effects and an evolution equation describing the spatiotemporal behavior 
of the convection is developed. It is found that in some parametric range, one convective 
cell forms across the layer like in the rigid-free case, while in other parametric domains, 
two convective cells emerge between the two free surfaces (which does not occur in the 
rigid-free case). 

The plan of the paper is as follows. Section II is devoted to the problem statement. 
Section III deals with a linear stability analysis of the phenomenon. In Section IV we 
derive the evolution equation describing the spatiotemporal behavior of the film undergoing 
Marangoni convection. Section V presents some further results of the study and discussion. 


II. PROBLEM STATEMENT 

A two-dimensional incompressible liquid layer of density p, viscosity p and thermal 
conductivity A is considered in this paper. It is assumed that this liquid layer is bounded 
vertically by two free surfaces with a thickness d between them and is of infinite extent in 
the horizontal direction. The ambient temperature of the air above the upper free surface 
and below the lower free surface is kept constant at T2 and Tj, respectively. We assume, 
without loss of generality since there is no gravity, that T2 < T\. The layer of interest, 
therefore, is exposed to a vertical temperature gradient —7 (7 > 0 ). Surface tension, a , 
acting at the interfaces is assumed to be temperature-dependent 
dai 


cr 1 — cr 1 0 + 


O’u — CT u0 + 


dT 

dcr u 

~dT 


{Tt-Tt 0 ) 
(T u - T u0 ) 
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where Ti o and T u o are the reference temperatures and the subscripts £, u correspond re- 
spectively to the lower and upper free surfaces of the layer. We note that, as is well-known, 
for most liquid-gas and liquid- liquid systems surface tension decreases with temperature, 
i.e. da j dT < 0, although this is not assumed in the following analysis. The ambient 
environment is considered to be passive. In what follows we focus on the case of a fluid 
layer with nondeformable static free boundaries. 

The governing equations describing the flow in the layer are: 

P(u t + UU X + VUy) — ~p x + fl(U XX + Uyy) 

p(v t + UV X + VVy) = -p y + p(v xx + Vyy) 

U X + Vy = 0 

Tf + uT X + vTy = k(T XX -f Tyy ) (l 0,6, C, <f) 

where ( u,v,p , T ) are respectively the velocity field components, pressure and temperature, 
and k is the thermal diffusivity of the fluid of interest. In what follows we consider the 
onset of Marangoni instability from the quiescent equilibrium state. Therefore it, v can be 
viewed as the velocity perturbations and T can be represented as 


T = T e - 7 (y - h t ) + 0(x,y,t) 


wherein is a perturbation of the temperature field. Eq.(ld) then can be rewritten 

as 


Ot + U0 X + Vdy - JV = k(9 XX + 0yy ) 


( 2 ) 


The boundary conditions at the lower, y - h t , and the upper, y = h u , free surfaces (both 
nondeformable) are: 
y = h t : 


(3a, b, c) 


pi^y ~\~ I'i) &l,x — 0 


X0 y — q e 0 = 0 


y = h u : 


v = 0 


p(u y + V x ) - a UiX = 0 
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Xdy + q u B — 0 6, c) 

Introducing d,d? j k, /c/d, ptc 2 /d 2 ,jd as scales for length, time, velocity, pressure and tem- 
perature, respectively, one rewrites eqs.(l)-(4) in the dimensionless form: 

U t + UU X + VUy = —p x + P(u xx + Uyy ) 


y = h t : 


V t + UV X + VVy = —Py + P(v xx + Vyy) 

Ux Vy 0 

B t +uB x + vB y-v = 0 xx + 0 y y (5 a,b,c,d) 


v = 0 


Uy + V x — MlO x = 0 


0 y -B t 0 = 0 


(6 ,a,b,c) 


V — h u : 


17 = 0 


Uy + V x + M U B x - 0 

By + B U B — 0 


(7a, 6, c) 


Here and from now on, the variables x,y,t,u,v,6 are respectively dimensionless spatial 
coordinates ( longitudinal and transverse), time, velocity perturbations (longitudinal and 
transverse) and temperature deviation from the conductive equilibrium state. 

The dimensionless parameters of the problem are: 

the Prandtl number P = 


the ’’lower” Marangoni number Mt = 
the ’’upper” Marangoni number M u = 


pi/K 

pl/K 


the ’’lower” Biot number Bt = 


the ’’upper” Biot number B u = 
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wherein qi,q u are the rates of heat transfer by convection at the lower and upper free 
surfaces, respectively. 

III. LINEAR STABILITY ANALYSIS 

In this section we study the stability limit of the quiescent state subject to infinitesimal 
perturbations, under the conditions that both free surfaces of the liquid layer are nearly 
insulating with regard to the temperature perturbations. An analytic expression of the 
critical values of the lower and upper Marangoni numbers at the onset of instability is 
derived. It is shown that, under certain conditions, the first unstable mode, i.e. the 
mode which is first amplified as the critical threshold is crossed, corresponds to a zero 
wavenumber, k=0. This behavior is similar to that of a rigid- free liquid layer with nearly 
insulating boundaries and indicates that the first instability is a longwave one. This fact 
will be used in the weakly nonlinear analysis in the next section. 

Our linear analysis follows in general the analysis made by Pearson 1 . Therefore, where 
possible, we will not go into full details. Linearizing the problem given by eqs. (5)-(7) one 
obtains: 

U% Px + ^(u xx d - ^yy) 

Ut Py d“ -P(^sex T v yy) 

U X + Vy = 0 

0 t —v = 0 XX + 0 yy (8a, 6, c, d ) 

with the boundary conditions 

at y — I 

v — 0 

Uy “ 1 “ v x “ 0 

0y = 0 (9a, fe, c) 

at y = 4 

v = 0 

u y + v x + M U 6 X — 0 
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(10a, b, c) 


e y = o 

Here we have assumed the boundaries of the fluid layer to be non-deformable, located at 
y — — | and y — ~2 an d nearly insulating (i.e. Bi,B u 0), which corresponds to the 
’’insulating case” considered by Pearson 1 . 

Eqs. (8a,b,c) can be reduced to 

v xxt + Vyyt = P{v X xxx + 2v xij ,3, + v yyyy ) (H) 

Introducing into eqs. (ll),(8d),(9),(10) 

v(x,y,t) = F(x)$(y)e ut 

0(x,y,t) = F(x)x(y)e wt (12) 

one obtains 

P(V" - 2k V' + fc 4 $) = u,($" - A: 2 *) (13) 

x" - (fc 2 + v)x - (1 4 ) 

and F" + k 2 F = 0, where k is the wavenumber of the perturbation in the x-direction. The 
boundary conditions corresponding to eqs. (9)-(10) are 

at y = — \ 


at V = -j 


$ = 0 , — k 2 Mix = 0 , x — 0 (15a, 6, c) 

$=0 , + k 2 M u x = 0 , x' = 0 (16a, £>, c) 


Looking for the conditions to be held at the neutral stability surface (u> = 0) one obtains 
the following solutions for eqs. (13), (14), (15a) and (16a) : 

2«S *5 

= ai [sinh(Aiy) y cosh(fcy)] + 02 [ cosh(fcy) + 2 y sinh(fcy)] (17) 

c c 

X = a 3 sinh(fcy) + cosh (ky) cosh(fcy) H — ^-(y cosh(/cy) — ky sinh(fcy))] + 

-y [-^ysinh(fcy) + -^(ysinh(fcy) - ky 2 cosh(fcy))] (18) 
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where s(k) _ sinh(A:/ 2 ) and c(fc) = cosh(A;/ 2 ). Eqs. (15b), (15c), (16c) are now used 
simultaneously with eqs. (17), (18) to determine the three unknowns a 1 ,a 2 ,a 3 . These are 
substituted into eq. (16b) to find the neutral stability criterion: 


where 


M u 


biMi + & 2 
b^Mi + b\ 


h = -8 csk* + 32 cVfc 2 
b 2 = 256c 3 s 3 k 4 

b3 = (c 2 + s 2 )k 3 - cs(k 4 + 3 A: 2 ) + 4 c 3 s 3 


(19) 


( 20 ) 


Figure 1 displays the behavior of M u as a function of Mi and k at criticality as 
expressed by eqs.(19)-(20). Typically in an experiment one would change both Marangoni 
numbers, such as by varying the temperature gradient 7 (see eq. (32)), the ratio of the 
lower and upper Marangoni numbers remaining nearly constant (since dcr/dT will remain 
nearly constant at both the top and bottom surfaces). Therefore it is useful to plot the 
critical curves M u — M u {k\a ), where 


Mi _ dai/dT 
Mu “ d^Jdf’ 


( 21 ) 


which result from intersections of the surface given by eq. (19) (see fig. 1 ) with the planes 
M t /M u = a for various values of a. Substituting M t = aM u in eq. (19) and solving for 
M u gives 


M« 


—bi(l - a) 


2b, 


3ft 


1 ± 4/1 + 


4&2&3CI 
52(1 - a )2 


( 22 ) 


where b u b 2 , and b 3 are given by eq. (20). (The functions 5i(Jb),5 2 (fc) and b 3 (k) are found 
to be positive and monotonically increasing as a function of k for k > 0 and vanish at 
k = 0 .) We note that for a = 0 , M u = b 2 /b 1 ] and that for a = 1 , M u = ± y /b^Jb 3 . 

Figure 2 shows some of these curves for various a. As can be seen from these curves, 
at onset the instability will be either to a finite wavenumber or to zero wavenumber. So 
an important question is: As M u and Mi are changed while keeping a constant, such as 
by increasing the temperature gradient 7 , under what conditions will the first instability 
be to a finite wavenumber and under what conditions will it be to zero wavenumber? To 
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find this condition, it is useful to look at the longwave limit (fc 
found to be (for - sign in eq. (22) and 1) 

24 


M u = 


8(3 — 7a + 3a 2 ) 
a 5(1 — a) 3 


0) of eq. (22) which is 


(23) 


The other root corresponding to the plus sign diverges as 360(1 — 1 ja)k~ 2 for k =* 0. Since 
we are mainly interested in the zero wavenumber instability, we will focus our attention 
on the root given by eq. (23). There are two primary cases to consider: a < 1 which at 
the onset of the longwave instability corresponds to M u > 0; and a > 1 which corresponds 
there to M u < 0. 

Case a < 1: If the coefficient of k 2 is negative in eq. (23), the critical curve M u (k;a) 
will initially decrease as k is increased from zero. It will then eventually increase as 
k is increased further, as given by eq. (22) and seen in fig. 2a. Therefore, as M u is 

gradually increased toward the critical curve from below, the first instability will be to 

finite wavenumber, since the minimum of the curve is at finite k. However, if the coefficient 
of k 2 is nonnegative, the critical curve M u (k,a) will increase as k is increased from zero. 
It then continues to increase monotonically as k is further increased as seen in figs. 2a, c. 
Therefore, as M u is gradually increased toward the critical curve from below, the first 
instability will be to zero wavenumber, since the minimum of the curve is at k = 0. 

Case a > 1: If the coefficient of k 2 is positive in eq. (23), the critical curve M u (k;a ) 
will initially increase as k is increased from zero. It will then eventually decrease as 
k is increased further, as given by eq. (22) and seen in fig. 2b. Therefore, as M u is 

gradually decreased toward the critical curve from above, the first instability will be to 

finite wavenumber, since the maximum of the curve is at finite k. However, if the coefficient 
of k 2 is nonpositive, the critical curve M u (k;a ) will decrease as k is increased from zero. 
It then continues to decrease monotonically as k is further increased as seen in fig. 2b. 
Therefore, as M u is gradually decreased toward the critical curve from above, the first 
instability will be to zero wavenumber, since the maximum of the curve is at k = 0. 

Examining the coefficient of k 2 in eq. (23) we thus find that the first instability will 
be to zero wavenumber if 


a < 1 — ^ 0.56574 or a > 7+ y^ = 1.76759 


in which case the onset of instability will occur at 


(24) 


M u -M t = 24 


(25) 


9 



where a is defined in eq. (21). Eq. (25) was found from taking fc = 0 in eq. (23) and 
replacing a by Mi/M u . ^From eq. (23) for k = 0 and eq. (24) we see that the limits on 
M u at a zero wavenumber instability are 

-31.2666 £ 12(1 - a/13) < M u < 12(1 + >/l3) £ 55.2666 (26) 

In contrast, the first instability will be to finite wavenumber when 

7 — \/l3 7 + a/13 

< a < — 

6 6 

We note that when dcr jdT is the same for the upper and lower surfaces, which corresponds 
to a — 1, M u diverges as 24\/l5 fc _1 as k =* 0 and has a minimum at finite k , corresponding 
to a finite wavenumber instability, 17 which is included in condition (27). We also note that 
the root corresponding to the plus sign in eq. (22), which as noted above diverges as k =± 0, 
corresponds to a finite wavenumber instability. 

In the limit k =± oo all neutral curves from eq. (22) are found to be given by 

M u ~ 8 k 2 or M u ~ — —k 2 (28) 

a 

the former giving the limiting behavior of the curves for which M u > 0 and the latter 
giving the limiting behavior of the curves for which M u < 0. The former is similar to the 
behavior found by Pearson for rigid-free boundary conditions. 1 

Figures 3a and 3b give the critical values of the Marangoni number and wavenumber, 
M tt)C and k c , respectively, at which instability first sets in as a function of a. These values 
correspond to the minimum or maximum (depending on whether M u > 0 or M u < 0) 
of the curves given in Fig. 2. For the case a < 0, as seen in Fig. 2c, only the lower of 
the two curves for a given value of a: is relevant, instability occurring when M u is above 
the minimum of the curve. M u being above the minimum of the upper of the two curves 
for a given a < 0 simply corresponds to there being a band of stable wavenumbers, the 
system as a whole still being unstable. The critical value of the upper Marangoni number, 
M u ,c exists for all values of a, when M u > 0. In contrast, Af u>c exists only for positive 
a increasing indefinitely for a =* 0 when M u < 0. Let us reiterate that as follows from 
our previous analysis and also from inspection of Fig. 3b, the first instability given by the 
M u > 0 branch is to zero wavenumber for ct < 0.56574 and to finite wavenumber otherwise. 
For the M u < 0 branch the first instability is to zero wavenumber for a > 1.76759 and to 
finite wavenumber otherwise. 
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A particular case of interest is the one when surface tension at the lower free surface 
is temperature independent: Mi — 0 (a = 0). Then the critical condition for the longwave 
instability is M u = 24. The corresponding critical Marangoni number for the case of the 
upper free surface and the lower rigid plane was found by Pearson 1 : M = 48. In our case 
(free-free) the critical Marangoni number is lower than in the Pearson’s (rigid-free) case. 
The reason is that a rigid boundary provides a certain damping for the emerging instability 
although generated at the opposite boundary. 


IV. WEAKLY NONLINEAR ANALYSIS. 


We turn now to the weakly nonlinear analysis of stability of a quiescent state of a 
fluid layer of uniform thickness (we neglect the deformability of the free surfaces in what 
follows) which is heated from below and is subject to surface tractions due to surface tension 
gradients. As stated before, we will focus on the case of nearly insulating free surfaces. 
The equilibri um temperature distribution corresponding to the quiescent conductive state 
in the fluid layer is 


T — T\ — B u (Ti - T 2 )~ 


Biy + 1 


(31) 

Bt + B u (l + B t ) y ’ 

Here we assumed the liquid layer to be contained between y = 0 and y = 1, for convenience. 
The equilibrium temperature gradient across the layer, 7 = (Ti — T u )/d , is therefore 


7 = 


B u B t {T x -T 2 ) 


d[Bi + I?u(l + Bt)\ 

Introduce the stretched spatial and temporal variables by 

C - X V* > V = V , t = e 2 i 


(32) 


(33) 


where e is a small perturbation parameter measuring the distance above onset [see eq. 
(36c, d)]. Note that the variable x is stretched, whereas the variable y is not stretched, 
which is consistent with a longwave instability slightly above onset [also see discussion 
following eq. (58)]. The Biot numbers of the system are assumed to be small and of order 


e 2 as 


B u = e 2 fd u , B t = e 2 /?/ 


(34) 
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The appropriate scalings for the streamfunction of the perturbation flow, 'f, and for the 
deviation of the temperature from the equilibrium, eq.(31) are 11 

ip = (35) 


0 = 0 

The scaled variables, then, and 6 are expanded in powers of e: 

V = i + e 2 V 2 + ... 


0 = 0o + e0i + e 2 02 + ... 

Mi — Mio + tMi\ 4- ... 

M u — M u o -j- tM u \ -J- ... (36a, 6, c, d) 

wherein Mio,M u q represent respectively the values of the lower and the upper Marangoni 
number at criticality. There is still an arbitrariness in the definition of e which can be 
removed by assuming a value for M u \ or Mn (or some linear combination of them). For 
example, taking M ul = M u0 gives e = (M u - M u0 )/M u0 . 

Eliminating the pressure terms from eqs.(5)-(7) and introducing the scaling (35) yields: 

e 3 tf CCr + e 2 V VT , T + - e^^ VVT , + = 


+ 2e*« vv + * vvvv) 

e 2 Qr + €$„0c - €$<0, + etf < = £0 a + 0„„ (37 a, b ) 


T] = 0 : 


V = 1 : 


= 0 


-€* <C ) + M*0 C =0 

0, - e 2 /?j0 = 0 


(38a, b, c ) 


V = 0 


(<S VV - etf <c ) + M u 0 c = 0 

0, + e 2 j3 u O = 0 


(39a, b, c) 
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Introducing the expansions ( 36 ) into eqs. ( 37 )-( 39 ) yields in the leading order 

^07777771; — 0 


©07777 = 0 

( 40 a, b) 

•CJ 

II 

O 

* 

O 

II 

O 


^0777; — Af, o 0 o < — 0 


©077 = 0 

( 41 a, b , c) 

*1 = 1 : 

^o=0 


^07777 + M u o©o^ = 0 


0 

O 

-=) 

II 

0 

( 42 a, b, c ) 

A solution of eqs. ( 40 b), ( 41 c), ( 42 c) is given by 


0 

0 

II 

( 43 ) 

where this function is as yet unknown. Further, an evolution equation in terms of /(£,r), 
describing the spatiotemporal behavior of the flow and the temperature fields, will be 
derived. A solution for the leading order of the streamfunction, 'Ifo? is given by 

*0 = (A z r ] 3 + A 2 tj 2 + A^)f 

( 44 ) 

where 

A\ = — ~Mio + —M u o , A 2 = —Mio , As — — — (M^ 0 
3 6 2 0 

+ M u o) 

and prime denotes differentiation with respect to £. 

The velocity field, therefore, at this order of approximation is 


uo = \£ot7 = (3A377 2 + 2A277 + Ai)f* 


Uo = -’J'oc = (A 3 tj 3 + A 2 T ) 2 + AiT])f" 

( 45 a, b) 

Along the free boundaries the velocities are 


0 

II 

0 

II 

0 

3 

0 

II 

pr- 

( 46 ) 
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u o — ( 3 A 3 4- 2 A 2 + A\ )f' , vq — 0 


77 = 1 : 


(47) 


Integrating eq.(37b) at the leading order between the free surfaces with respect to 77 and 
using eqs. (43), (44) one obtains the condition to be held at criticality: 


M u 0 — Mi 0 = 24 


(48) 


which coincides with the relation found using linear stability analysis in the previous 
section, eq.(25). 

At the first order of approximation one obtains 


^ lr ) T ) i rn vv 


©17777 — — © 0 « + ^077©0< + ^OC 


(49a, b) 


with the boundary conditions 
77 = 0 : 


^ 1=0 

— ^17777 + -^<i©oc + Mm©i< + ^occ = 0 

01 , =0 


(50a, b, c) 


77 = 1 : 


^ 1=0 

^17777 + M u i ® o <; + M u o®i{ — = 0 

©i, = 0 

The solution for the temperature is given by 


(51a, 6 , c) 


©i = -^/' 2 (6 tW + 4-4, r, 3 + 3-4,7,-*)+ -i/"(-30, 2 + 10-4, ,, 3 + 5-4 2 7j 4 + 34.,, 5 ) + g((, r) 

(52) 

where q(C,r) is an arbitrary function. The streamfunction is determined via 

\Tr, — ~^ 3 f n M _L I ( ^2 fifii ^3 - 

Wl " 70P ;/ V + 30 P f} 71 +( lo p ff )v + 

( 22 p g “/ ,,, ) r ? 4 + 93V 3 + (Mnf + M io q ')-y + g ^77 (53) 
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where 


f 1 Q* 

g\ — (M u i — 2 Mn)— + (M u o — 2Mto) — + 

, M u o A-\ j4_2 7A| 2t4. 2 ^3 ^^-3 w/ r» . 

V o ' mn ' nn r> ' i r D ' OK T} J* ^ 1 


3 ' 12P 90P 

4j 4. 2 2.4.3 


15P 


35P 


1 S/f / 1 . -^1 I 4 2 , -^3 ^ rill 

— + u o( Y 2 + ~36" + 72 + 120 


( 54 ) 


and 

53 


Ma + M u i r , Mio + M u o , ( M u 0 j4i4 2 j4 2 , 4 2 A3 43 ^ , fU 

= * — / « — <? -(-T- + -^- + ^ + "^“ + : r7i^^ } 


6 ' 6 3 v 3 6P ' 9P 6P 10P 


r A-i 

+1_ 6" + 2 


4 2 4.3 M u 0 ,, 4i 


+ 4^ + 


; (i- 


•4 2 43 . , ni 


■)\f 


(55) 


2 ' 12 v 3 6 10 

To simplify eqs. (54), (55); eq. (48) and the ensuing identity A.i/2-t- 4. 2 /3 + A$/4 = 1 was 
used. 

The evolution equation for the temperature disturbance f(C, T ) is derived by integrat- 
ing eq.( 37 b) at the order e 2 with respect to 77 across the layer: 

fr+[ — ^0<Ql77 + , ® r l< — ®icc) + (fiu + fil)f — 0 (56) 

Jo 

Introducing eqs.(44), (52), (53) into eq.(56) yields the key result of the work 
Sr + "if"" + 7T2 /" + fSf ~ 7 T 3 (/' 3 )' + 7 T 4 (/' 2 )" = 0 


(57) 


where the coefficients are given by 


f3 — fit + /?u 


7Tl 


+ 


7T 2 = 


M u0 M 2 0 
360 8640 

M u i - Mn 


*3 = 


105 


(128 


24 

M, 


(58) 


uo , M 2 0 


+ 


72 


) 


11 . 1 

Tl = T0 ( ~2 + 




7 P' v 12 

It is readily shown that the coefficient 7 Ti is positive in the range 12(1 — \/l3) < M u o < 
12(1 + v / 13) which constitutes the range of validity of the present nonlinear theory and 
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which coincides with eq.(26) for the range of M u for which onset is to zero wavenumber. 
Moreover, 7 r 3 is always positive. On the other hand, the 7r 2 -term provides instability to 
the considered system, if 7r 2 > 0, i.e. M ui - M tl > 0. In this case, the difference between 
the upper and the lower Marangoni numbers, AM, is 


AM — M u — Mi — M u o — Mio + e(M a i — Mu ) > 24 


and the layer of interest is unstable with respect to longwave perturbations, in the specified 
range of M u o as also follows from the linear stability analysis brought in the previous 
section. 

Also, as follows from eqs. (57), (58) and independently from the linear analysis, the 
linear stage of the evolution of an infinitesimal temperature disturbance is not affected by 
the value of the Prandtl number, P, which appears only in the nonlinear 7 r 4 - term. 

Eq. (57) and its various particular cases have recently appeared in the literature 
s, n,i4, i »-22 R e 04 is particularly devoted to a study of eq. (57). By rescaling the spatial, 
temporal and dependent variables, £,r, / respectively, 


/■ > 1/ 2 - 1/2 _ 2 , , 1/2 - 1/2 
C =* C^l ^2 > T =± T7Ti7r 2 , f f IT/ 7 T 3 ' 


one obtains the equation 


fr + /"" + f" + sf- (f«y + u = 0 (59) 

in the domain 0 < < < L wherein 

s = /?7Ti7r^“ 2 , u; = 2tt±'K 1 
We note that eq. (59) can be recast into the form 


fr + [(1 - f n )}‘\' + /"" + Sf + u(/' 2 )" = 0 (60) 


The cubic term / ,3 in eq. (59) is shown to ensure the emergence of bounded amplitude 
patterns. Moreover, in the case of uj — 0 corresponding to M u q = 12, eq. (59) is symmetric 
under the transformation / — / and has a Lyapunov functional bounded from below (if 

•s > 0) such that eq. (59) can be rewritten in the form: 


fr = 


6F 

Tf 


(61) 
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with a free energy functional given by 


r -fji* , + 1 S-±r + ±rv' 

and A being either the periodic or infinite domain. 


(62) 


V. FURTHER RESULTS AND DISCUSSION 

Turning back to the leading order of the streamfunction \foj eq.(44), one finds that 
the polynomial .A377 3 + A 2 T ] 2 + A\T] has three roots: 


2M (0 — M u 0 M u 0 — 48 

1)1 - 0 , 712 - ,m- M>0 + M(0 - 2JW u o — 24 


(63) 


For M u 0 > 48 or M u 0 < —24, 773 lies in the domain 0 < 77 < 1, and therefore, two 
different cells across the layer are expected to form, one below the other. This can be 
also understood from the fact, that for M u 0 > 48 both the velocities at the free surfaces, 
uq(t) = 1) and uo(rj — 0), given by eqs. (46), (47) are negative when /' > 0, and positive 
otherwise. For M u0 < —24 both u 0 (t? = 0) and uo(tj = 1) are positive for positive /' > 0, 
and negative otherwise. 

In contrast, for —24 < M u 0 < 48 the velocities at the free boundaries have different 
signs, leading to formation of a single convective cell across the layer between the adjacent 
extrema of /(£,t). The effect of the emergence of two convective cells across the layer is 
novel and is not observed in the ’’rigid- free” case 8 . 

It also readily follows from eq.(63) that for M u 0 > 48 the lower cell is smaller than the 
upper one. If M u 0 < —24 the upper cell is smaller that the lower one. These considerations 
are observed in our numerical study of eq.(57) as well. 

Equation (57) is solved numerically using a time-splitting method together with the 
periodic boundary conditions in the domain 0 < C < L = 57 r. The value of the Prandtl 
number is chosen as P = 5 which corresponds to water. The distance from the criticality 
expressed by the difference M u i — Mu is taken as 6. Equation (57) is amended with the 
initial condition related to the fundamental mode /o(C) = 0.01 sin(27r£/L). 

Figure 4a presents the solution for eq.(57), /, the leading order of the temperature 
disturbance, for M u o — 54,M^o == 30, (3 = 0. The corresponding flow field is displayed in 
Fig. 4b. [The streamlines shown in Figs. 4b, 5b, and 6b were calculated using eq. (44).] As 
explained above, two convective cells emerge across the layer, the lower cell being smaller 
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than the upper one. The latter can be understood from the fact that the strength of 
surface tractions at the upper surface is higher than that at the lower one, \M u q | > \Miq\. 
The direction of the flow at both surfaces (see eqs. (46,47)) is from the hottest spot at the 
surfaces corresponding to the maximum of / (see fig. 4a) to the coldest one corresponding 
to the minimum of /, in agreement with the direction of the surface shear stress arising 
from the local variation of surface tension. Note that the temperature of the fluid is given to 
the leading order by T(y) + /(£, r), where the first term is the linear equilibrium state, (eq. 
(31)) and the second one is obtained from solving eq. (57). At this order of approximation, 
therefore, the temperature varies coherently along the free boundaries and the hottest and 
the coldest spots on each of them are located at the same vertical cross-section. 

Figure 5a shows the solution for eq. (57) with M u0 = —30 , M^o = — 54, /? = 0. The 
corresponding flow field is displayed in Fig. 5b. In this case, \Miq\ > \M U o|, and the lower 
cell is larger than the upper one. The direction of the flow at both surfaces is from the 
coldest spot at the surfaces corresponding to the minimum of / toward the hottest one 
corresponding to the maximum of / (see fig. 5a). 

A third type of flow is presented in Figs. 6a and 6b. The parameter values taken here 
are Mio = — 12, M u q = 12 , (3 = 0. This set of parameters corresponds to 7c± = 0, therefore 
the solution obtained is symmetric with respect to the midplane. The steady pattern for 
the temperature disturbance is displayed in the first and the flow field which consists of 
one convective cell across the layer is displayed in the second. The direction of the flow 
along the lower free surface is from the coldest spot to the hottest one (due to the negative 
Marangoni number) and the flow along the upper free surface is from the hottest spot to 
the coldest one. 

In this paper we have studied the onset of the Marangoni instability in a liquid layer 
vertically bounded by two free nearly insulating surfaces with arbitrary interfacial prop- 
erties. Using the linear stability analysis, we found the stability ( instability) domains 
in the parametric space including the disturbance wavenumber and the both Marangoni 
numbers. In particular we derived the criterion for the first instability being to zero or 
to finite wavenumber. The critical threshold for the onset of the longwave instability is 
found to be given by the relationship between the upper and lower Marangoni numbers 
: M u — Mi ~ 24. A weakly nonlinear analysis is applied to study the spatiotemporal 
evolution of the flow and temperature fields at the onset of the longwave instability. It 
was shown that in some parametric range one convective cell forms across the layer while 
in other parametric domains two convective cells emerge between the two free surfaces. 
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APPENDIX: THE CASE OF ARBITRARY BIOT NUMBERS B t ,B u 

When the Biot numbers at the free surfaces, B t and B u , are not both small the critical 
threshold for the instability is given by [compare eq. (19)] 

6i Mi + b 2 

Af u — 7 Tj 1 r* 

b$Mi + 64 


where 


*>! = -8 csk 4 - 8 B u {c 2 + s 2 )k 3 + 8 cs(B u + 4 c 2 s 2 )k 2 + 16 B u c 2 s 2 {c 2 + s 2 )k 

b 2 = 256cVife 4 + 128 (B e + B u )c 2 s 2 (c 2 + s 2 )k 3 + 256 B t B u c 3 s z k 2 
b 3 = (c 2 + s 2 )k 3 - cs{k 4 + 3k 2 ) + 4 c 3 s 3 
b 4 = —Scsk 4 - 8 Bt(c 2 + s 2 )k 3 + 8 cs(B t + 4c 2 s 2 )A: 2 + 16R^c 2 ^ 2 (c 2 + s 2 )k 

Again solving for M u (k',a.) and expanding about k — 0 we find that there is a root 
that is finite at k = 0 only if B t + B u + B t B u = 0 or if B e = B u = 0. Therefore, there is no 
longwave instability for any combination of Biot numbers and a weakly nonlinear stability 
analysis of the type derived in this paper is not relevant (except for B t =* 0 and B u =* 0, 
as studied in the main body of the paper). 
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Figure 1 .—The neutral surface M u - k) as given by 

eq. (19). In the case of longwave perturbations, k => 0, the 
cross-section of the surface is the straight line given by 
M u — Mf = 24, eq. (25). The instability domain is above 
the surface. 



k 



k 


Figure 2.— The neutral stability curves for various values of 
the parameter a, as given by eq. (22). Two curve 
branches correspond to each value of a ^ 0. (a) 0 < a < 1 . 

If M u < 0 the first instability is to finite wavenumber 
perturbations. For M u > 0 the first instability is to zero 
wavenumber perturbations when a < 0.56574 and to finite 
wavenumber otherwise, (b) a > 1 . For M u > 0 the first 
instability is to finite wavenumber perturbations. For M u < 0 
the first instability is to zero wavenumber perturbations for 
a > 1 .76759 and to finite wavenumber otherwise, (c) a. <> 0. 

The first instability is the longwave one. Note that the unstable 
domain lies between the two curves corresponding to the 
given value of a. There is only one curve corresponding to 
a = 0. 
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Figure 3. — (a) The critical value of the upper Marangoni 
number, M u c , as a function of the parameter a. Note 
that in the case of the negative Marangoni numbers the 
corresponding curve with the opposite sign is displayed, 
(b) The critical wavenumber of the perturbation, fc c , at 
which the first instability occurs. Note that in the case 
of M u > 0, k c = 0 for a < 0.56574, while k c = 0 for oc> 
1.76759 when M u < 0. 



Figure 4. — The steady-state solution: (a) the temperature 
deviation from the equilibrium state; (b) the streamlines 
of the fluid flow, for eq. (57) with periodic boundary 
conditions, the initial condition fo(£) = 0-^ s\n{2'nt/L) > 
and P = 5, p = 0, - M n = 6, L = 5tt, = 54, 

M/o = 30. The direction of the flow in the cells at the free 
surfaces is from the hottest spot to the coldest one. 

The streamlines in the lower cell are not displayed due 
to its thinness. The two convective cells are counter- 
rotating. 







Figure 5.— The steady-state solution: (a) the temperature 
deviation from the equilibrium state and (b) the streamlines 
of the fluid flow, for eq. (57) with periodic boundary 
conditions, the initial condition = 0.01 sin(2ir C/L), and 
P = 5,P = 0,M^-M n =6,{. = Sir, Myo = -30, M yo = -54. 
The direction of the flow in the cells at the free surfaces is 
from the coldest spot to the hottest one. The streamlines 
in the upper cell are not displayed due to its thinness. 

The two convective cells are counterrotating. 



Figure 6.— The steady-state solution: (a) the temperature 
deviation from the equilibrium state; (b) the streamlines 
of the fluid flow, for eq. (57) with periodic boundary cond- 
itions, the initial condition f 0 (Q = 0.01 stn(2irJyL), and P = 5, 
0 = 0, - M n = 6, L = 5, Myo = 12, = -12. The 

direction of the flow in the cells at the upper free surface is 
from the hottest spot to the coldest one, and at the lower 
free surface from the coldest spot to the hottest one. The 
streamlines of the flow are symmetric with respect to the 
midplane of the layer. 
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